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Abstract We consider the computation of Bernoulli, Tangent (zag), and Secant (zig 
or Euler) numbers. In particular, we give asymptotically fast algorithms for comput- 
ing the first n such numbers in 0(n 2 (\ogn) 2+ °^) bit-operations. We also give very 
short in-place algorithms for computing the first n Tangent or Secant numbers in 
0(n 2 ) integer operations. These algorithms are extremely simple, and fast for mod- 
erate values of n. They are faster and use less space than the algorithms of Atkinson 
(for Tangent and Secant numbers) and Akiyama and Tanigawa (for Bernoulli num- 
bers). 



1 Introduction 

The Bernoulli numbers are rational numbers B„ defined by the generating function 

(1) 



n>0 



n\ exp(z) — 1 



Bernoulli numbers are of interest in number theory and are related to special val- 
ues of the Riemann zeta function (see jj2). They also occur as coefficients in the 
Euler-Maclaurin formula, so are relevant to high-precision computation of special 
functions Q §4.5]. 
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It is sometimes convenient to consider scaled Bernoulli numbers 

(2n)! 

with generating function 

LC„z 2 " = -^-. (3) 
tanh(z/2) 

The generating functions ((TJ and (O only differ by the single term B\z, since the 
other odd terms vanish. 

The Tangent numbers T n , and Secant numbers S„, are defined by 

2n-l 2n 

T, T "n TT7 =taiu ' E 5 "7TTT = secz - < 4 ) 

In this paper, which is based on an a talk given by the first author at a workshop 
held to mark Jonathan Borwein's sixtieth birthday, we consider some algorithms for 
computing Bernoulli, Tangent and Secant numbers. For background, combinatorial 
interpretations, and references, see Abramowitz and Stegun Q] Ch. 23] (where the 
notation differs from ours, e.g. (— l) n E2 n is used for our S„), and Sloane's |29l 
sequences A000367, A000182, A000364. 

Let M(n) be the number of bit-operations required for 77-bit integer multipli- 
cation. The Schonhage-Strassen algorithm 1271 gives M(n) = (9(nlognloglog«), 
and Fiirer ifTTl has recently given an improved bound M(n) = (9(n(logn)2 log "). 
For simplicity we merely assume that M(n) = 0(n(\ogn) 1+0 ^>), where the o(l) 
term depends on the precise algorithm used for multiplication. For example, if 
the Schonhage-Strassen algorithm is used, then the o(l) term can be replaced by 
log log log n / log log « . 

In §*2r[3]we mention some relevant and generally well-known facts concerning 
Bernoulli, Tangent and Secant numbers. 

Recently, Harvey |20l showed that the single number B n can be computed in 
0(« 2 (log«) 2+o W) bit-operations using a modular algorithm. In this paper we show 
that all the Bernoulli numbers Bq, . . . , B„ can be computed with the same complexity 
bound (and similarly for Secant and Tangent numbers). 

In |2] we give a relatively simple algorithm that achieves the slightly weaker 
bound 0(« 2 (log«) 3+o ( 1 '). In fJ3]we describe the improvement to <9(« 2 (logn) 2+ °' 1 '). 
The idea is similar to that espoused by Steel |3T1 , although we reduce the problem to 
division rather than multiplication. It is an open question whether the single number 
Bin can be computed in o(n 2 ) bit-operations. 

In S|6]we give very short in-place algorithms for computing the first n Secant or 
Tangent numbers using 0(n ) integer operations. These algorithms are extremely 
simple, and fast for moderate values of n (say n < 1000), although asymptotically 
not as fast as the algorithms given in §[0]42] Bernoulli numbers can easily be de- 
duced from the corresponding Tangent numbers using the relation (fl4] i below. 
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2 Bernoulli Numbers 



From the generating function (Q]) it is easy to see that the B„ are rational numbers, 
with Z?2n+i = if n > 0. The first few nonzero B n are: Bq = l,B\ = — 1/2, Z?2 = 1/6, 
fl 4 = -1/30, B 6 = 1 /42, fig = -1/30, Bio = 5/66, Bi 2 = -691/2730, B H =7/6. 

The denominators of the Bernoulli numbers are given by the Von Staudt - 
Clausen Theorem lfT2"l[3"0"l . which states that 

B' 2n :=B 2n + £ \ 

[p-\)\2n P 

Here the sum is over all primes p for which p — 1 divides 2n. 

Since the "correction" B' 2n — B 2 „ is easy to compute, it might be convenient in a 
program to store the integers B' 2n instead of the rational numbers B 2n or C„. 

Euler found that the Riemann zeta-function for even non-negative integer argu- 
ments can be expressed in terms of Bernoulli numbers - the relation is 



N „_i B 2 „ = 2C(2«) 
(2ft)! (2tt) 2 



Since £(2n) = 1 + 0(4 ") as n -> +oo, we see that 

2(2n)! 



|B 2n | 



(27T 



i2« 



From Stirling's approximation to (2ft)!, the number of bits in the integer part of 
B 2n is 2ft lg n + 0(n) (we write lg for log 2 ). Thus, it takes £2 (n 2 log n) space to store 
B\,...,B n . We can not expect any algorithm to compute B\,...,B n in fewer than 
£2(n 2 \ogn) bit-operations. 

Another connection between the Bernoulli numbers and the Riemann zeta- 
function is the identity 

£*--«-> (6 > 

for n £ Z, « > 1. This follows from (O and the functional equation for the zeta- 
function, or directly from a contour integral representation of the zeta-function 11331 . 

From the generating function ([T), multiplying both sides by exp(z) — 1 and equat- 
ing coefficients of z, we obtain the recurrence 

k /k+l\ 

£ \ Bj = for k>0. (7) 

This recurrence has traditionally been used to compute Bq, . . . ,B 2n with 0(n 2 ) arith- 
metic operations, for example in 11221 . However, this is unsatisfactory if floating- 
point numbers are used, because the recurrence is numerically unstable: the relative 
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error in the computed Z?2« is of order A"e if the floating-point arithmetic has preci- 
sion e, i.e. lg(l/e) bits. 

Let C n be defined by (0. Then, multiplying each side of by sinh(z/2)/(z/2) 
and equating coefficients gives the recurrence 

X Cj i 

(8) 



H {2k +1-2 j) ! 4 k ~J {2k) ! A k 

Using this recurrence to evaluate Cq,C\ ,C n , the relative error in the computed 
C„ is only 0(n 2 e), which is satisfactory from a numerical point of view. 

Equation (0 can be used in several ways to compute Bernoulli numbers. If we 
want just one Bernoulli number Bi n then £(2n) on the right-hand-side of (0 can be 
evaluated to sufficient accuracy using the Euler product: this is the "zeta-function" 
algorithm for computing Bernoulli numbers mentioned (with several references to 
earlier work) by Harvey 1201 . On the other hand, if we want several Bernoulli num- 
bers, then we can use the generating function 

^ = -2 £ (-l)W*, (9) 

computing the coefficients of z 2k , k < n, to sufficient accuracy, as mentioned in 
l9l [l0l . This is similar to the fast algorithm that we describe in jJU The similarity can 
be seen more clearly if we replace nz by z in (0, giving 

_i^ = - 2 y(-ifM z a ( io) 

tanh(z) ^to n 2k 

since it is the rational number £(2«)/7z; 2 " that we need in order to compute B2,, 
from 0. In fact, it is easy to see that ( TUH i is equivalent to (0. 

There is a vast literature on Bernoulli, Tangent and Secant numbers. For exam- 
ple, the bibliography of Dilcher and Slavutskii lfT31 contains more than 2000 items. 
Thus, we do not attempt to give a complete list of references to related work. How- 
ever, we briefly mention the problem of computing irregular primes I0QT], which 
are odd primes p such that p divides the class number of the p-th cyclotomic field. 
The algorithms that we present in §Sj4]-0below are not suitable for this task because 
they take too much memory. It is much more space-efficient to use a modular algo- 
rithm where the computations are performed modulo a single prime (or maybe the 
product of a small number of primes), as in (8][TTl[T4j|20l. Space can also be saved 
by the technique of "multisectioning", which is described by Crandall fPH §3.2] and 
Hare Q9). 
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The Tangent numbers T„ (n > 0) (also called Zag numbers) are defined by 

In— 1 

L z z " . sinz 

fI> o r "(^ =tanZ = ^' 

Similarly, the Secant numbers S„ (n > 0) (also called Euler or Zig numbers) are 
defined by 

z 2 " 1 
£ S n — — = secz = 

n% 0) ! C °SZ 

Unlike the Bernoulli numbers, the Tangent and Secant numbers are positive integers. 
Because tanz and secz have poles at z = 7t/2, we expect T„ to grow roughly like 
(2« — 1 ) ! (2 /n)" and S„ like (2n) ! (2 / n) n . To obtain more precise estimates, let 

CbW = (l-2-')C(*) = l+3-' + 5-* + - 
be the odd zeta-function. Then 

T„ 2 2 "+ l ^(2n) 2 2 " +1 



(2n-l) 



(11) 



(this can be proved in the same way as Euler's relation (01 for the Bernoulli num- 
bers). We also have 0Q (23.2.22)] 

S n _2 2 "+ 2 j3(2« + l) 2 2 "+ 2 



(2n)\ K 2n+l n 2n+l 

where 

P(s) = i{-l) J (V d3) 

7=0 

From (0 and (fTTT i, we see that 

T„ = (-1)"- 1 2 2 "(2 2 "-1)— ■ (14) 
2n 

This can also be proved directly, without involving the zeta-function, by using the 
identity 

1 2 

tanz = — — ■ 

tanz tan(2z) 

Since T n G Z, it follows from ( fPfl i that the odd primes in the denominator of Bi„ 
must divide 2 2 " — 1 . This is compatible with the Von Staudt-Clausen theorem, since 
(p — l)\2n implies p|(2 2 " — 1) by Fermat's little theorem. 

T„ has about An more bits than \B2„] , but both have 2n lg n + 0(n) bits, so asymp- 
totically there is not much difference between the sizes of T„ and |~B2«1 ■ Thus, if our 



6 



Richard P. Brent and David Harvey 



aim is to compute B^„, we do not lose much by first computing T n , and this may be 
more convenient since T„ € Z, Bin £ Q. 



4 A Fast Algorithm for Bernoulli Numbers 

Harvey l20l showed how B„ could be computed exactly, using a modular algorithm 
and the Chinese remainder theorem, in (9(« 2 (logn) 2+0 0) bit-operations. The same 
complexity can be obtained using (0 and the Euler product for the zeta-function 
(see the discussion in Harvey lEUl §1]). 

In this section we show how to compute all of Bq, . . . ,B„ with almost the same 
complexity bound (only larger by a factor (9(log«)). In §3]we give an even faster 
algorithm, which avoids the <9(log«) factor. 

Let A(z) = «o + o.\z + aiz 1 H be a power series with coefficients in R, with 

«o 7^ 0. Let B(z) = bo + biz-\ be the reciprocal power series, so A(z)B(z) = 1. 

Using the FFT, we can multiply polynomials of degree n— 1 with O(nlogn) real 
operations. Using Newton's method [24 28 1, we can compute bo,... ,b„-i with the 
same complexity 0(n\ogn), up to a constant factor. 

Taking A(z) = (exp(z) — l)/z and working with A^-bit floating-point numbers, 
where = nlg(n) + 0(n), we get Bq,...,B„ to sufficient accuracy to deduce 
the exact (rational) result. (Alternatively, use (|3} to avoid computing the terms 
with odd subscripts, since these vanish except for B\.) The work involved is 
(9(«log«) floating-point operations, each of which can be done with A^-bit ac- 
curacy in 0(n(logn) 2+o W) bit-operations. Thus, overall we get Bq,...,B„ with 
(9(« 2 (log«) 3+ "(^) bit-operations. Similarly for Secant and Tangent numbers. We 
omit a precise specification of and a detailed error analysis of the algorithm, since 
it is improved in the following section. 



5 A Faster Algorithm for Tangent and Bernoulli Numbers 

To improve the algorithm of 2] for Bernoulli numbers, we use the "Kronecker- 
Schonhage trick" ||7] §1.9]. Instead of working with power series A(z) (or poly- 
nomials, which can be regarded as truncated power series), we work with binary 
numbers A (z) where z is a suitable (negative) power of 2. 

The idea is to compute a single real number si which is defined in such a way 
that the numbers that we want to compute are encoded in the binary representation 
of si . For example, consider the series 



k>Q 
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The right-hand side is an easily-computed rational function of z, say A(z). We use 
decimal rather than binary for expository purposes. With z = 10~ 3 we easily find 

1001000 

A 10~ 3 ) = = 0.001004009016025036049064081 100 •■ • 

v ' 997002999 ~ ~ 

Thus, if we are interested in the finite sequence of squares (1 2 ,2 2 ,3 2 , . . . , 10 2 ), it is 
sufficient to compute srf = A(10~ 3 ) correctly rounded to 30 decimal places, and we 
can then "read off" the squares from the decimal representation of srf . 

Of course, this example is purely for illustrative purposes, because it is easy to 
compute the sequence of squares directly. However, we use the same idea to com- 
pute Tangent numbers. Suppose we want the first n Tangent numbers (T\ , T2, . . . , T n ). 
The generating function 



tanz = £ T k 

k>\ 



z 2k ~ l 
(2k-\)\ 



gives us almost what we need, but not quite, because the coefficients are rationals, 
not integers. Instead, consider 

(2fl-l)!tam = £ T[ n z 2k - 1 +R n (z), (15) 

k=\ 



where 



(2«-l)! 
(2*-l)! 



TL= "—, /.'. Tk (16) 



is an integer for 1 < k < n, and 



2k-\ 

Rn(z)= Z K n z 2k - l = (2n-l)\ £ T kl ±—- (17) 

k=n+\ k=n+\ \ LK l l- 



is a remainder term which is small if z is sufficiently small. Thus, choosing 
z = 2~ p with p sufficiently large, the first 2np binary places of (2n — l)!tanz de- 
fine T( n , T 2 ' n , . . . , T„ . Once we have computed T[ n , n , . . . , T„ n it is easy to deduce 
7i,r 2 ,'. . . ,f„ from 

ji 

T k,n 
Ik 



(2n-l)!/(2fe-l)! 
For this idea to work, two conditions must be satisfied. First, we need 

0< T[ n < l/z 2 = 2 2p , 1 <k<n, (18) 

so we can read off the Tl from the binary representation of (2« — l)!tanz. Since 
we have a good asymptotic estimate for 7J., it is not hard to choose p sufficiently 
large for this condition to hold. 

Second, we need the remainder term R„(z) to be sufficiently small that it does 
not influence the estimation of T' „. A sufficient condition is 
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< R„(z) < z 2 "- 1 . (19) 

Choosing z sufficiently small (i.e. p sufficiently large) guarantees that condition ( fT9l 
holds, since R„(z) is 0(z 2n+i ) as z — >• with n fixed. 

Lemmas[3]and|4]below give sufficient conditions for (fT8T i and (fT9l to hold. 

Lemma 1. 

/ 2 n2(H) 

< - for/t>l. 



(2*-l) 
Proof. From (JTTJ, 

7* / 2 \ ^ / 2 ^ ^ / 9 \ ^ 7r2 / 9 \ 1) 



2 - 5,(2*) <2 ~ Co(2)< - - = - □ 



(2*-l)! 

Lemma 2. (2« - 1)! < n Zn ~ l for n > 1. 

(2b - 1)! =n n(n - ;)(« + ;) = » ft* 2 - f) < n 2 "- 1 
7=1 7=1 

with equality iff « = 1 . □ 

Lemma 3. Ifk >l,n>2,p = \n\g(n)\, z = 2~ p , and T' kn is as in (fT6l l. then z < n~" 
ant/ r/ n < 1 /z 2 . 

Proof. We have z = 2~P = 2- r«ig(»)l < 2-» 1 s(») = n ~", which proves the first part 
of the Lemma. 

Assume k > 1 and n > 2. From LemmaQ] we have 



^<(2»-l)l(|) 
and from Lemma|2]it follows that 



2(*-l) 

<(2n-l)\, 



U n <n 2 "- y <n 2 ". 

From the first part of the Lemma, n 2 " < l/z 2 , so the second part follows. □ 

Lemma 4. If n > 2, p = [~nlg(n)~|, z = 2~ p , and R„(z) is as defined in ( 1171 i. then 
0<R n {z) <0.lz 2 "- 1 . 

Proof. Since all the terms in the sum defining R„ (z) are positive, it is immediate that 
R„{z)>0. 
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Since « > 2, we have p > 2 and z < 1/4. Now, using LemmaQ] 



k=n+l 



<(2n-l)! E - 

k=n+l \' L 



2(H) 



2k-l 



\ 2 / n \ 4 

2 \ 2 «+i / , . / 2z\ . /2z 



<(2n-l)!(-] [!+(-= 



In I / \ 2 



2 A 2„+i / / . 1 2z 



<(2n-l)!^-] z 

Since z < 1/4, we have 1/(1 - (2z/7r) 2 ) < 1 .026. Also, from LemmaEl 
(2n — 1)! < n 2 " -1 . Thus, we have 

^ (Z) < 1.026„ 2 "-^ 2 ^" 



z 2„-l 



Now z 2 < n 1,1 from the first part of Lemma [3] so 

The right-hand side is a monotonic decreasing function of n, so is bounded above 
by its value when n = 2, giving Rn(z)/Z 2 "- 1 <0.1. □ 

A high-level description of the resulting Algorithm FastTangentNumbers is given 
in Figure Q] The algorithm computes the Tangent numbers 7\ , Tj , . . . , T„ using the 
Kronecker-Schonhage trick as described above, and deduces the Bernoulli numbers 
#2,^4, • • • ,B2n from the relation ( TBI . 



Input: integer n > 2 

Output: Tangent numbers T\,...,T n and (optional) Bernoulli numbers 62,-84, • • - ,B2„ 

p*- r«igwi 

S^Eo<t<„(-l)^ +1 x (2n)!/(2*+l)! 
C^Z 0ik<n (-l) k z 2k x(2n)\/(2k)\ 

V <— [z l ~ 2n x (2« — 1) ! x S/C] (here |*] means round x to nearest integer) 

Extract T' kn = T^(2« — l)\/(2k— 1)!, 1 <k<n, from the binary representation of V 

T k ^T^ n x(2k-iy./(2n-l)\,k = n,n-l y ...,l 

Bik <- ('- 1 )*" 1 (A x 71 /2 2 ^ 1 ) / (2 2k -\),k=\,2,...,n (optional) 

return T\ , T2, ■ ■ ■ , T n and (optional) B2, B4, . . . , B2n 



Fig. 1 Algorithm FastTangentNumbers (also optionally computes Bernoulli numbers) 
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In order to achieve the best complexity, the algorithm must be implemented 
carefully using binary arithmetic. The computations of S (an approximation to 
(2n) ! sinz) and C (an approximation to (2ft) ! cosz) involve computing ratios of fac- 
torials such as (2»)!/(2fe)!, where < k < n. This can be done in time (9(« 2 (log«) 2 ) 
by a straightforward algorithm. The iV-bit division to compute S/C (an approxima- 
tion to tanz) can be done in time 0(Nlog(N) loglog(A^)) by the Schonhage-Strassen 
algorithm combined with Newton's method §4.2.2]. Here it is sufficient to take 
N = 2np + 2 = 2rc 2 lg(«) + 0(n). Note that 

y = £ 2 2(n-k)p T ^ (21) 
k=l 

is just the finite sum in (H5]> scaled by z 1 " 2 " (a power of two), and the integers 
can simply be "read off" from the binary representation of V in n blocks of 
2/7 consecutive bits. The T' kn can then be scaled by ratios of factorials in time 

0{n 2 (\ogn) 2+ "^) to give the Tangent numbers 7i,72,... , T„. 

The correctness of the computed Tangent numbers follows from Lemmas [3}j4] 
apart from possible errors introduced by S/C being only an approximation to tan(z). 
Lemma[5]shows that this error is sufficiently small. 

Lemma 5. Suppose that n > 2, z, S and C as in Algorithm FastTangentNumbers. 
Then 



z 1 - 2 "(2«-l)! 



S 

tanz 

C 



< 0.02. (22) 



Proof. We use the inequality 



B' 



< 



|A|-|B-B'| + |B|-|A-A'| 



\B\-\B' 



(23) 



Take A = sinz, B = cosz, A' = S/{2n)\, B' = C/(2ft)! in (|23]>. Since n > 2 we have 
< z < 1/4. Then |A| = |sinz| < z. Also, |S| = cosz > 31/32 from the Taylor 
series cosz = 1 — z 2 /2 + ■ ■ ■ , which has terms of alternating sign and decreasing 
magnitude. By similar arguments, \B'\ > 31/32, \B— B'\ < z 2 " /(2n)!, and 
| A — A' | < z 2 " +1 /(2« + 1)!. Combining these inequalities and using (1231) . we obtain 



S 

--tanz 



6-32-32z 2 " +I 1.28z 2 " +1 

< - ~7Z — TT < 



5-31-31 (2ft)! 



(2n)! 



Multiplying both sides by z 1 2 "(2« — 1 ) ! and using 1.28z 2 /(2«) <0.02, we obtain 
the inequality (1221) . This completes the proof of Lemma[5] □ 

In view of the constant 0.02 in d22l and the constant 0. 1 in Lemma|H the effect of 
all sources of error in computing z'~ 2 "(2n — 1)! tanz is at most 0.12 < 1/2, which is 
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too small to change the computed integer V, that is to say, the computed V is indeed 
given by d2TT i. 

The computation of the Bernoulli numbers #2,^4, • • • >^2n from T\,...,T n , is 
straightforward (details depending on exactly how rational numbers are to be repre- 
sented). The entire computation takes time 

O(jV(logA0 1+o(1) ) = 0(n 2 (logn) 2+o «). 
Thus, we have proved: 

Theorem 1. The Tangent numbers T\,...,T n and Bernoulli numbers Bi,B^, . . . ,B% n 
can be computed in 0(n 2 (logn) 2+0 0) bit-operations using 0(n 2 \ogn) space. 

A small modification of the above can be used to compute the Secant numbers 
So, Si,. . . ,S„ in <9(« 2 (log«) 2+0 '^) bit-operations and (9(« 2 log«) space. The bound 
on Tangent numbers given by Lemma[T]can be replaced by the bound 

(2«)! - 

which follows from (TTSl i since j3 (2« + 1 ) < 1 . 

We remark that an efficient implementation of Algorithm FastTangentNumbers 
in a high-level language such as Sage ll32l or Magma Q is nontrivial, because 
it requires access to the internal binary representation of high-precision integers. 
Everything can be done using (implicitly scaled) integer arithmetic - there is no 
need for floating-point - but for the sake of clarity we did not include the scaling 
in Figure Q] If floating-point arithmetic is used, a precision of N bits is sufficient, 
where N = 2np + 2. 

Comparing our Algorithm FastTangentNumbers with Harvey's modular algo- 
rithm ll20l . we see that there is a space-time trade-off: Harvey's algorithm uses less 
space (by a factor of order n) to compute a single B„, but more time (again by a fac- 
tor of order n) to compute all of B\, . . . ,B n . Harvey's algorithm has better locality 
and is readily parallelisable. 

In the following section we give much simpler algorithms which are fast enough 
for most practical purposes, and are based on three-term recurrence relations. 



6 Algorithms Based On Three- Term Recurrences 

Akiyama and Tanigawa 121] gave an algorithm for computing Bernoulli numbers 
based on a three-term recurrence. However, it is only useful for exact computations, 
since it is numerically unstable if applied using floating-point arithmetic. It is faster 
to use a stable recurrence for computing Tangent numbers, and then deduce the 
Bernoulli numbers from (TT~4b - 
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6.1 Bernoulli and Tangent numbers 



We now give a stable three-term recurrence and corresponding in-place algorithm 
for computing Tangent numbers. The algorithm is perfectly stable since all opera- 
tions are on positive integers and there is no cancellation. Also, it involves less arith- 
metic than the Akiyama-Tanigawa algorithm. This is partly because the operations 
are on integers rather than rationals, and partly because there are fewer operations 
since we take advantage of zeros. 

Bernoulli numbers can be computed using Algorithm TangentNumbers and the 
relation dT4b . The time required for the application of ( TBI is negligible. 

The recurrence (l24l that we use was given by Buckholtz and Knuth |23l , 
but they did not give our in-place Algorithm TangentNumbers explicitly. Related 
recurrences with applications to parallel computation were considered by Hare (19\ . 

Write t = tan*, D = d/dx, so Dt = 1 +t 2 and D(t") = nt n ~ l {l+t 2 ) for n > 1. 
It is clear that D"t is a polynomial in t, say P n {t). Write P n {t) = T,j>oPnjt^- Then 
deg(P„) = n + 1 and, from the formula for D(t n ), 

Pn.j = (;'- l)p„-ij_i + (;' + l)p„_ij+i. (24) 

We are interested in 7^ = (d/dx) 2k ~ l tanx |. v= n = /?2t_i(0) = P2fc-i,0> which can be 
computed from the recurrence (l24l i in 0(k 2 ) operations using the obvious boundary 
conditions. We save work by noticing that p n j = if n + j is even. The resulting 
algorithm is given in Figure |2] 

Input: positive integer n 

Output: Tangent numbers T\,...,T n 

for k from 2 to n 

T k ^(k-l)T k ^ 
for k from 2 to n 
for j from k to n 

7)«-C/-*)rv-i + (/-* + 2)2) 
return T\,T2, . ■ ., T„, 



Fig. 2 Algorithm TangentNumbers 

The first for loop initializes 3j = Pk-\,k = — !)!■ The variable is then used 
to store pk.k-i, Pk+ik-2> ■ ■ •> p2k-2 i> P2fc-i o at successive iterations of the second 
for loop. Thus, when the algorithm terminates, 7). = P2fc-i,o> as expected. 

The process in the case « = 3 is illustrated in Figure [5] where 7l' m ^ denotes the 
value of the variable 7^ at successive iterations m = 1,2,... ,«. It is instructive to 
compare a similar Figure for the Akiyama-Tanigawa algorithm in lETl . 

Algorithm TangentNumbers takes @(« 2 ) operations on positive integers. The in- 
tegers T n have (9(«log«) bits, other integers have (9(log«) bits. Thus, the overall 
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Ti =Po,i 
/ \ 

„(i) 

Ti = Pi.o = pi i2 

\ / \ 

T (2) _ „ T (l) _ „ 

s \ s 

T (2) _ n T (2) _ „ 

'2 — Pifi ^3 — P3,2 

\ / 
t(3) 

/ 

^3 = P5.0 

Fig. 3 Dataflow in Algorithm TangentNumbers for n = 3 

complexity is (9(n 3 (log«) 1+0 (^) bit-operations, or (9(« 3 logn) word-operations if n 
fits in a single word. 

The algorithm is not optimal, but it is good in practice for moderate values of n, 
and much simpler than asymptotically faster algorithms such as those described in 
§fflS For example, using a straightforward Magma implementation of Algorithm 
TangentNumbers, we computed the first 1000 Tangent numbers in 1.50 sec on a 
2.26 GHz Intel Core 2 Duo. For comparison, it takes 1.92 sec for a single A^-bit 
division computing T in Algorithm FastTangentNumbers (where Af = 19931568 
corresponds to n = 1000). Thus, we expect the crossover point where Algorithm 
FastTangentNumbers actually becomes faster to be slightly larger than n = 1000 
(but dependent on implementation details). 



6.2 Secant numbers 



A similar algorithm may be used to compute Secant numbers. Let s = sec x, t = tan.Y, 
and D = d/dx. Then Ds = st, D 2 s = s(l +2f 2 ), and in general D"s = sQ„(t), where 
Q„(t ) is a polynomial of degree n in t . The Secant numbers are given by Sk = Q2k(Q)- 
Let Q„(t) = Lt>o From 

D{st k )=st k+l +kst k -\\+t 2 ) 

we obtain the three-term recurrence 

q n +i,k = kq„.k-\ + (k + 1 )q ni k+ 1 for \<k<n. (25) 

By avoiding the computation of terms q n # that are known to be zero (« + k odd), and 
ordering the computation in a manner analogous to that used for Algorithm Tangent- 
Numbers, we obtain Algorithm SecantNumbers (see Figure|4|i, which computes the 
Secant numbers in place using non-negative integer arithmetic. 
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Input: positive integer n 

Output: Secant numbers So, S\ . . . , S n 

for k from 1 to n 

Sk <- fcS*_i 
for k from 1 to n 

for j from k + 1 to n 

Sj<-(J-k)Sj-i + U-k+i)Sj 
return Sq,Si , . . . ,S„. 



Fig. 4 Algorithm SecantNumbers 



6.3 Comparison with Atkinson's algorithm 

Atkinson J2] gave an elegant algorithm for computing both the Tangent numbers 
7*i , ?2, . . . , T„ and the Secant numbers Sq,S\,...,S„ using a "Pascal's triangle" style 
of algorithm that only involves additions of non-negative integers. Since a trian- 
gle with In + 1 rows in involved, Atkinson's algorithm requires 2n 2 + 0(n) integer 
additions. This can be compared with n 2 /2 + 0(n) additions and « 2 + 0{n) multi- 
plications (by small integers) for our Algorithm TangentNumbers, and similarly for 
Algorithm SecantNumbers. 

Thus, we might expect Atkinson's algorithm to be slower than Algorithm Tan- 
gentNumbers. Computational experiments confirm this. With n = 1000, Algorithm 
TangentNumbers programmed in Magma takes 1.50 seconds on a 2.26 GHz Intel 
Core 2 Duo, algorithm SecantNumbers also takes 1.50 seconds, and Atkinson's al- 
gorithm takes 4.51 seconds. Thus, even if both Tangent and Secant numbers are re- 
quired, Atkinson's algorithm is slightly slower. It also requires about twice as much 
memory. 
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